/**
 * @file exponentials.h
 * @brief Math approximations for computing exponentials.
 * @date January 16, 2019
 * @author Guillaume Giudicelli, MIT, Course 22 (g_giud@mit.edu)
 */

#ifndef EXPONENTIALS_H_
#define EXPONENTIALS_H_

/**
 * @brief Computes the F1 exponential term.
 * @details This method computes (1-exp(-x))/x using a 7th order polynomial
 *          approximation. It is accurate to 1e-5 over [0, 5]. It was generated
 *          by Colin Josey using Remez's algorithm.
 * @param x the value at which to evaluate the function, usually tau
 * @param expF1 a pointer to the exponential
 */
inline void expF1_poly(FP_PRECISION x, FP_PRECISION* expF1) {

  const FP_PRECISION p0 = 1.;
  const FP_PRECISION p1 = -4.998618823537523 * 1E-1;
  const FP_PRECISION p2 = 1.660264339632089 * 1E-1;
  const FP_PRECISION p3 = -4.0607961247327 * 1E-2;
  const FP_PRECISION p4 = 7.459558151235148 * 1E-3;
  const FP_PRECISION p5 = -9.773063624328772 * 1E-4;
  const FP_PRECISION p6 = 8.004982165323072 * 1E-5;
  const FP_PRECISION p7 = -3.016010858852142 * 1E-6;

  *expF1 = p7*x + p6;
  *expF1 = *expF1*x + p5;
  *expF1 = *expF1*x + p4;
  *expF1 = *expF1*x + p3;
  *expF1 = *expF1*x + p2;
  *expF1 = *expF1*x + p1;
  *expF1 = *expF1*x + p0;
}


/**
 * @brief Computes the F2 exponential term.
 * @details This method computes (x-2+exp(-x)*(2+x))/x**2 using a 8th order
 *          polynomial approximation. It is accurate to 1e-5 over [0, 5]. It was
 *          generated by Colin Josey using Remez's algorithm.
 * @param x the value at which to evaluate the function, usually tau
 * @param expF2 a pointer to the exponential
 */
inline void expF2_poly(FP_PRECISION x, FP_PRECISION* expF2) {

  const FP_PRECISION p1 = 1.666656766985825 * 1E-1;
  const FP_PRECISION p2 = -8.331262914972137 * 1E-2;
  const FP_PRECISION p3 =  2.492325839109710 * 1E-2;
  const FP_PRECISION p4 = -5.440953156443790 * 1E-3;
  const FP_PRECISION p5 =  9.034802244154581 * 1E-4;
  const FP_PRECISION p6 = -1.091608727341864 * 1E-4;
  const FP_PRECISION p7 =  8.411465095972204 * 1E-6;
  const FP_PRECISION p8 = -3.029020287833148 * 1E-7;

  *expF2 = p8*x + p7;
  *expF2 = *expF2*x + p6;
  *expF2 = *expF2*x + p5;
  *expF2 = *expF2*x + p4;
  *expF2 = *expF2*x + p3;
  *expF2 = *expF2*x + p2;
  *expF2 = *expF2*x + p1;
  *expF2 = *expF2*x;
}


/**
 * @brief Computes the H exponential term.
 * @details This method computes (1-exp(-x)*(1+x))/x**2 using a 8th order
 *          polynomial approximation. It is accurate to 1e-5 over [0, 5]. It was
 *          generated by Colin Josey using Remez's algorithm.
 * @param x the value at which to evaluate the function, usually tau
 * @param expH a pointer to the exponential
 */
inline void expH_poly(FP_PRECISION x, FP_PRECISION* expH) {

  const FP_PRECISION p0 = 0.5;
  const FP_PRECISION p1 = -3.33307480097059  * 1E-1;
  const FP_PRECISION p2 =  1.248597190564637 * 1E-1;
  const FP_PRECISION p3 = -3.305711771115656 * 1E-2;
  const FP_PRECISION p4 =  6.667701492411682 * 1E-3;
  const FP_PRECISION p5 = -1.028726890908856 * 1E-3;
  const FP_PRECISION p6 =  1.145870114989106 * 1E-4;
  const FP_PRECISION p7 = -8.079532805720403 * 1E-6;
  const FP_PRECISION p8 =  2.657707693467421 * 1E-7;

  *expH = p8*x + p7;
  *expH = *expH*x + p6;
  *expH = *expH*x + p5;
  *expH = *expH*x + p4;
  *expH = *expH*x + p3;
  *expH = *expH*x + p2;
  *expH = *expH*x + p1;
  *expH = *expH*x + p0;
}


/**
 * @brief Computes an exponential G term, used to reconstruct F1, F2 and H.
 * @details This method computes 1/x-(1-exp(-x))/x**2 using a 5/6th order
 *          rational approximation. It is accurate to 2e-7 over [0, 1e5].
 *          However, accuracy when reconstructing F1, F2 and H is lower.
 *          It was generated by Colin Josey using Remez's algorithm.
 * @param x the value at which to evaluate the function, usually tau
 * @param expG a pointer to the exponential
 */
inline void expG_fractional(FP_PRECISION x, FP_PRECISION* expG) {

  /* Coefficients for numerator */
  const FP_PRECISION p0 = 0.5;
  const FP_PRECISION p1 = 1.76558112351595 * 1E-1;
  const FP_PRECISION p2 = 4.041584305811143 * 1E-2;
  const FP_PRECISION p3 = 6.178333902037397 * 1E-3;
  const FP_PRECISION p4 = 6.429894635552992 * 1E-4;
  const FP_PRECISION p5 = 6.064409107557148 * 1E-5;

  /* Coefficients for denominator */
  const FP_PRECISION d0 = 1.0;
  const FP_PRECISION d1 = 6.864462055546078 * 1E-1;
  const FP_PRECISION d2 = 2.263358514260129 * 1E-1;
  const FP_PRECISION d3 = 4.721469893686252 * 1E-2;
  const FP_PRECISION d4 = 6.883236664917246 * 1E-3;
  const FP_PRECISION d5 = 7.036272419147752 * 1E-4;
  const FP_PRECISION d6 = 6.064409107557148 * 1E-5;

  FP_PRECISION num, den;
  den = d6*x + d5;
  den = den*x + d4;
  den = den*x + d3;
  den = den*x + d2;
  den = den*x + d1;
  den = den*x + d0;
  den = 1.f/den;

  num = p5*x + p4;
  num = num*x + p3;
  num = num*x + p2;
  num = num*x + p1;
  num = num*x + p0;

  *expG = num*den;
}


/**
 * @brief Computes the F1 exponential term.
 * @details This method computes (1-exp(-x))/x using a 5/6th order
 *          rational approximation. It is accurate to 2e-6 over [0, 1e5].
 *          It was generated by Colin Josey using Remez's algorithm.
 * @param x the value at which to evaluate the function, usually tau
 * @param expF1 a pointer to the exponential
 */
inline void expF1_fractional(FP_PRECISION x, FP_PRECISION* expF1) {

  /* Coefficients for numerator */
  const FP_PRECISION p0 = 1.0;
  const FP_PRECISION p1 = 2.4172687328033081 * 1E-1;
  const FP_PRECISION p2 = 6.2804790965268531 * 1E-2;
  const FP_PRECISION p3 = 1.0567595009016521 * 1E-2;
  const FP_PRECISION p4 = 1.0059468082903561 * 1E-3;
  const FP_PRECISION p5 = 1.9309063097411041 * 1E-4;

  /* Coefficients for denominator */
  const FP_PRECISION d0 = 1.0;
  const FP_PRECISION d1 = 7.4169266112320541 * 1E-1;
  const FP_PRECISION d2 = 2.6722515319494311 * 1E-1;
  const FP_PRECISION d3 = 6.1643725066901411 * 1E-2;
  const FP_PRECISION d4 = 1.0590759992367811 * 1E-2;
  const FP_PRECISION d5 = 1.0057980007137651 * 1E-3;
  const FP_PRECISION d6 = 1.9309063097411041 * 1E-4;

  FP_PRECISION num, den;

  den = d6*x + d5;
  den = den*x + d4;
  den = den*x + d3;
  den = den*x + d2;
  den = den*x + d1;
  den = den*x + d0;
  den = 1.f / den;

  num = p5*x + p4;
  num = num*x + p3;
  num = num*x + p2;
  num = num*x + p1;
  num = num*x + p0;

  *expF1 = num*den;
}


/**
 * @brief Computes the F2 exponential term.
 * @details This method computes (x-2+exp(-x)*(2+x))/x**2 using a 5/6th order
 *          rational approximation. It is accurate to 2e-6 over [0, 1e5].
 *          It was generated by Colin Josey using Remez's algorithm.
 * @param x the value at which to evaluate the function, usually tau
 * @param expF2 a pointer to the exponential
 */
inline void expF2_fractional(FP_PRECISION x, FP_PRECISION* expF2) {

  /* Coefficients for numerator */
  const FP_PRECISION p1 = 1.666661470036759 * 1E-1;
  const FP_PRECISION p2 = 3.59041632356318 * 1E-2;
  const FP_PRECISION p3 = 7.675127136944033 * 1E-3;
  const FP_PRECISION p4 = 6.408491755085618 * 1E-4;
  const FP_PRECISION p5 = 1.367575707015872 * 1E-4;

  /* Coefficients for denominator */
  const FP_PRECISION d0 = 1.0;
  const FP_PRECISION d1 = 7.153333128932897 * 1E-1;
  const FP_PRECISION d2 = 2.541555663123697 * 1E-1;
  const FP_PRECISION d3 = 5.613392571426973 * 1E-2;
  const FP_PRECISION d4 = 9.476002327852898 * 1E-3;
  const FP_PRECISION d5 = 9.145637477815584 * 1E-4;
  const FP_PRECISION d6 = 1.367575707015872 * 1E-4;

  FP_PRECISION num, den;
  num = p5*x + p4;
  num = num*x + p3;
  num = num*x + p2;
  num = num*x + p1;
  num = num*x;

  den = d6*x + d5;
  den = den*x + d4;
  den = den*x + d3;
  den = den*x + d2;
  den = den*x + d1;
  den = den*x + d0;
  *expF2 = num/den;
}


/**
 * @brief Computes the H exponential term.
 * @details This method computes (1-exp(-x)*(1+x))/x**2 using a 5/7th order
 *          rational approximation. It is accurate to 2e-6 over [0, 1e5].
 *          It was generated by Colin Josey using Remez's algorithm.
 * @param x the value at which to evaluate the function, usually tau
 * @param expH a pointer to the exponential
 */
inline void expH_fractional(FP_PRECISION x, FP_PRECISION* expH) {

  /* Coefficients for numerator */
  const FP_PRECISION p0 = 0.5;
  const FP_PRECISION p1 = 5.599412483229184 * 1E-2;
  const FP_PRECISION p2 = 1.294939509305754 * 1E-2;
  const FP_PRECISION p3 = 2.341166644220405 * 1E-3;
  const FP_PRECISION p4 = 3.686858969421769 * 1E-5;
  const FP_PRECISION p5 = 4.220477028150503 * 1E-5;

  /* Coefficients for denominator */
  const FP_PRECISION d0 = 1.0;
  const FP_PRECISION d1 = 7.787274561075199 * 1E-1;
  const FP_PRECISION d2 = 2.945145030273455 * 1E-1;
  const FP_PRECISION d3 = 7.440380752801196 * 1E-2;
  const FP_PRECISION d4 = 1.220791761275212 * 1E-2;
  const FP_PRECISION d5 = 2.354181374425252 * 1E-3;
  const FP_PRECISION d6 = 3.679462493221416 * 1E-5;
  const FP_PRECISION d7 = 4.220477028150503 * 1E-5;

  FP_PRECISION num, den;
  num = p5*x + p4;
  num = num*x + p3;
  num = num*x + p2;
  num = num*x + p1;
  num = num*x + p0;

  den = d7*x + d6;
  den = den*x + d5;
  den = den*x + d4;
  den = den*x + d3;
  den = den*x + d2;
  den = den*x + d1;
  den = den*x + d0;
  *expH = num/den;
}


/**
 * @brief Computes the G2 exponential term.
 * @details This method computes 2/3 - (1 + 2/x) * (1/x + 0.5 - (1 + 1/x) *
                                 (1-exp(-x)) / x) using a 5/5th order
 *          rational approximation. It is accurate to 1e-6 over [0, 1e6].
 *          It was generated by Colin Josey using Remez's algorithm.
 * @param x the value at which to evaluate the function, usually tau
 * @param expG2 a pointer to the exponential
 */
inline void expG2_fractional(FP_PRECISION x, FP_PRECISION* expG2) {

  /* Coefficients for numerator */
  const FP_PRECISION a1 = -8.335775885589858 * 1E-2;
  const FP_PRECISION a2 = -3.603942303847604 * 1E-3;
  const FP_PRECISION a3 = 3.7673183263550827 * 1E-3;
  const FP_PRECISION a4 = 1.124183494990467 * 1E-5;
  const FP_PRECISION a5 = 1.6837426505799449 * 1E-4;

  /* Coefficients for denominator */
  const FP_PRECISION b1 = 7.454048371823628 * 1E-1;
  const FP_PRECISION b2 = 2.3794300531408347 * 1E-1;
  const FP_PRECISION b3= 5.367250964303789 * 1E-2;
  const FP_PRECISION b4 = 6.125197988351906 * 1E-3;
  const FP_PRECISION b5 = 1.0102514456857377 * 1E-3;

  FP_PRECISION num, den;
  num = a5 * x + a4;
  num = num * x + a3;
  num = num * x + a2;
  num = num * x + a1;
  num *= x;

  den = b5 * x + b4;
  den = den * x + b3;
  den = den * x + b2;
  den = den * x + b1;
  den = den * x + 1.f;

  *expG2 = num / den;
}


/**
 * @brief Evaluates the function 1 - exp(x) for x negative. Vectorises well.
 *        Accurate to 6.18 digits (single precision) for entire domain (including
 *        near zero, unlike intrinsic)
 *
 *        Valid for x ~= [-1.5e6, 0]
 * @param x value at which to evaluate 1 - exp(x)
 * @param expv value of function
 */
inline void cram7(FP_PRECISION x, FP_PRECISION* expv) {

  /* Generated in Mathematica, accurate to 6.18 digits (single precision),
     tau [-1.5e6,0] */
  const FP_PRECISION c1n = -1.00000014302666667201396424463;
  const FP_PRECISION c2n = 2.34841040052684510704433796447 * 1E-1;
  const FP_PRECISION c3n = -6.24785939603762121316592924635 * 1E-2;
  const FP_PRECISION c4n = 1.00434102711342948752684759736 * 1E-2;
  const FP_PRECISION c5n = -1.35724435934263932676353754751 * 1E-3;
  const FP_PRECISION c6n = 9.51474224366003625378414851577 * 1E-5;
  const FP_PRECISION c7n = -1.60076055315534285575266516209 * 1E-5;

  const FP_PRECISION c0d = 1.;
  const FP_PRECISION c1d = -7.34847118148952339633322706422 * 1E-1;
  const FP_PRECISION c2d = 2.63193362386411901729092564316 * 1E-1;
  const FP_PRECISION c3d = -6.09467155163113059870970359654 * 1E-2;
  const FP_PRECISION c4d = 1.00863490579686697359577926719 * 1E-2;
  const FP_PRECISION c5d = -1.35667018708833025497446407598 * 1E-3;
  const FP_PRECISION c6d = 9.51502816434275317085698085885 * 1E-5;
  const FP_PRECISION c7d = -1.60076032420105715765981718742 * 1E-5;

  FP_PRECISION num, den;

  den = c7d;
  den = den * x + c6d;
  den = den * x + c5d;
  den = den * x + c4d;
  den = den * x + c3d;
  den = den * x + c2d;
  den = den * x + c1d;
  den = den * x + c0d;

  num = c7n;
  num = num * x + c6n;
  num = num * x + c5n;
  num = num * x + c4n;
  num = num * x + c3n;
  num = num * x + c2n;
  num = num * x + c1n;
  num = num * x;

  *expv = num / den;
}


#endif // EXPONENTIALS_H_
